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Preface 


This  report  describes  the  research  carried  out  during  six  month  of  2008-2009 
in  the  Robotics  and  Manufacturing  Automation  (RAMA)  Laboratory  of  Duke 
University.  The  major  effort  during  this  time  was  the  development  of  a  theoret¬ 
ical  approach  for  using  samples  of  stochastic  processes  for  solving  the  stochastic 
optimal  control  problems.  This  effort  is  described  in  the  following  sections. 

During  this  time  period  a  technical  paper  dealing  with  force  field  estimation 
based  on  video  data  was  also  prepared  and  submitted  for  publication  considera¬ 
tion  in  a  Special  Issue  on  Physical  System  Modeling  of  the  Journal  of  Dynamic 
Systems,  Measurement,  and  Control,  Series  G  of  the  ASME  Transactions.  The 
submitted  paper’s  abstract  is  attached  in  Appendix  A.  Publication  decision 
on  this  paper  is  still  awaited.  A  scale-down  copy  of  the  poster  presented  at 
the  CMPI  Symposium  on  Multi-Scale  Modeling  of  Host/Pathogen  Interactions, 
June  23  -  25,  2009,  Pittsburgh,  PA,  is  attached  as  Appendix  B. 

The  financial  support  provided  by  the  United  States  Army  Research  Office 
for  the  Short  Term  Innovative  Research  (STIR)  project  entitled  ’’Stochastic 
Model-Based  Control  of  Multi  Robot  Systems”  under  Grant  Number  W911NF- 
08-1-0503  is  gratefully  acknowledged. 


1 


Abstract 

In  this  report,  we  consider  control  of  single-  and  multi-robot  systems  as  an 
optimal  control  problem.  Solution  of  this  problem  may  be  of  enormous  com¬ 
plexity  because  of  a  large-number  of  robots,  a  large  number  of  redundant  states, 
and  environmental  uncertainties.  Motivated  by  estimation  methods  based  on 
statistical  sampling  and  employed  for  solving  complex  estimation  problems,  we 
explore  the  possibility  of  using  stochastic  process  samples  for  computing  the 
optimal  control. 

In  our  work,  the  state  of  an  individual  robot  is  described  by  the  state  com¬ 
posed  of  discrete  and  continuous  state  variables.  We  model  the  robot  using  a 
hybrid  automaton  and  corresponding  system  of  partial  differential  equations  de¬ 
scribing  the  robot’s  state  probability  density  function.  For  solving  the  optimal 
control  problem,  we  use  the  minimum-principle  for  partial  differential  equations 
to  obtain  the  Hamiltonian  that  characterizes  the  optimal  control.  However,  hav¬ 
ing  in  mind  that  the  Hamiltonian  evaluation  based  on  the  partial  differential 
equations  is  computationally  expensive,  we  propose  and  explore  strategies  in 
which  the  Hamiltonian  evaluation  is  based  on  samples  from  stochastic  processes 
generated  from  the  individual  robot  model.  We  also  show  that  if  the  analysis  is 
limited  to  the  so-called  process  self-adjoint  dynamical  systems,  the  Hamiltonian 
evaluation  can  be  simplified. 

Using  computational  statistical  methods  opens  the  possibility  to  solve  con¬ 
trol  problems  in  robotics  in  real-time  within  the  stochastic  optimal  control 
framework.  This  possibility  also  depends  on  processor  design  that  implements 
computational  methods.  While  the  design  of  such  processors  is  a  separate  issue, 
ideally  it  should  be  driven  by  stochastic  processes,  such  as  noises  in  electronic 
components.  This  concept  can  ultimately  provide  miniature  computational 
hardware  for  solving  complex  multi-robot  control  problems,  in  which  compu¬ 
tations  are  driven  by  laws  of  statistical  physics. 
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-  set  of  real  numbers,  set  of  real  number  vectors  of  dimension  n. 

-  set  of  discrete  states,  i.e.,  integer  indexes  {1,  2, 3, . . .} 

-  set  of  admissible  control. 

-  probability  density  function  (PDF)  of  the  hybrid  state  at  time  t.  This 
variable  is  a  vector  of  functions,  it  depends  on  x  G  X  and  t  G  R,  but  x 
and  t  are  frequently  omitted  in  expressions. 

-  the  PDF  component  corresponding  to  the  discrete  state  i,  i  G  Q. 

-  the  adjoint  state  distribution 

-  the  discrete  approximation  of  the  adjoint  state  distribution. 

-  the  adjoint  state  PDF 

-  the  probability  of  the  discrete  state  i,  i  G  Q. 

-  the  probability  of  the  discrete  adjoint  state  i,  i  G  Q. 

-  the  transition  rate  matrix. 

-  the  component  of  the  transition  rate  matrix  [Ft]ij  =  Xij. 

-  the  transition  rate  matrix  that  depends  on  control  vector  u,  u  G  Uad- 

-  the  component  of  the  linear  operator  F  corresponding  to  the  vector  fields 
fi  of  discrete  states  i  G  Q. 

-  the  PDF,  the  control  and  time  dependent  Hamiltonian,  p,  u  and  t  are 
frequently  omitted. 

-  the  optimal  control 

-  the  expected  value  with  respect  to  the  state  PDF  p. 
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Stochastic  Model- Based  Control  of 


Multi-Robot  Systems 


1  Introduction 

A  multi-robot  system  is  a  collection  of  robots  designed  to  perform  a  specified 
task.  Collection  of  all  relevant  variables  uniquely  specifying  the  individual  robot 
state,  such  as  position,  status,  actuator  powers  etc.,  we  call  the  state.  In  this 
respect,  a  solution  of  a  multi-robot  control  problem  assumes  that  the  update 
of  each  individual  robot  state  is  in  accordance  with  the  specified  multi-robot 
system  task.  The  update  is  defined  by  the  specific  control  approach  or  ar¬ 
chitecture.  In  this  research,  we  have  developed  an  approach  in  which  the  task 
specification  is  provided  in  the  form  of  a  cost  function  that  should  be  maximized 
(or  minimized)  by  the  multi-robot  system.  Such  cost  functions  can  relate  to  the 
presence  of  robots  at  specific  locations,  or  to  the  amount  of  information  that 
can  be  gathered.  It  may  also  include  other  factors  such  as  energy  consumption 
or  amount  of  communicated  information. 

From  this  perspective,  the  multi-robot  system  control  can  be  considered  as 
an  optimal  control  problem  resulting  in  the  update  of  individual  robot  states 
leading  towards  a  minimum  (or  maximum)  of  the  cost  function  defining  the 
multi-robot  mission  task.  Solution  of  this  problem  may  be  of  enormous  com¬ 
plexity  because  of  a  large  number  of  robots,  a  large  number  of  redundant  states. 
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as  well  as  existing  environment  uncertainties.  In  this  report,  we  discuss  an 
approach  to  deal  with  this  complexity  computationally  for  both  single-  and 
multiple-robot  systems. 

For  many  years  it  has  been  known  that  the  optimal  control  and  optimal 
estimation  problems  are  dual  [6].  For  example,  we  use  the  optimal  control 
theory  to  derive  linear  quadratic  regulator  (LQR),  and  in  the  same  theoretical 
framework  we  can  derive  the  Kalman  filter  (KF).  In  the  first  case  we  design  the 
controller  that  minimizes  the  quadratic  cost  function  over  the  state  space,  and  in 
the  second  case  we  design  the  estimator  minimizing  the  cost  function  involving 
state  space  estimations.  Having  this  in  mind,  it  is  expected  that  estimation 
methods  based  on  statistical  sampling,  employed  for  solving  complex  estimation 
problems  [12],  can  be  also  exploited  for  solving  complex  control  problems  for 
single  and  more  important  multi-robot  systems  under  presence  of  uncertainty. 

Along  this  idea,  Kappen  et.  al.  [7,  13]  used  stochastic  differential  equations 
to  model  individual  agents.  Based  on  this  description,  it  is  possible  to  re¬ 
late  Hamilton- Jacob i-Bellman  partial  differential  equation  with  samples  of  the 
stochastic  process  trajectories  and  use  the  samples  to  define  the  stochastic  op¬ 
timal  control  of  multi-agent  systems.  In  this  framework,  the  state  of  individual 
robots  is  continuous.  However,  the  state  of  real  robots  is  generally  described 
by  a  combination  of  continuous  and  discrete  variables,  i.e.,  by  a  hybrid  state. 
Therefore,  it  is  more  natural  to  describe  the  robot  behavior  using  a  hybrid  au¬ 
tomaton  [2].  The  automaton  describes  the  discrete  and  continuous  variables 
change  in  time,  which  depends  on  events  influencing  the  robot  behavior.  In 
the  case  of  a  large-size  multi-robot  population,  it  becomes  highly  complex  to 
predict  events  from  the  robot  local  environment.  Because  of  that,  we  model  a 
large-size  robots’  population  considering  the  stochastic  hybrid  model  and  study 
how  it  can  be  controlled. 
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In  this  research  project,  we  consider  a  problem  in  which  the  presence  of  a 
large-size  robot  population,  in  a  desired  region  of  operating  space,  is  maximized. 
This  problem  is  formulated  in  the  hybrid  system  framework  in  [11].  Its  solution, 
based  on  the  minimum  principle  for  partial  differential  equations,  is  presented 
in  [9,  10],  but  it  is  solved  numerically  only  in  the  case  when  the  presence  of 
the  robots  is  maximized  along  one  dimension  (ID).  The  difficulty  in  developing 
more  realistic  2D  control  is  rooted  in  the  fact  that  that  the  Hamiltonian,  which 
defines  the  optimal  control,  includes  integral  terms  that  depend  on  the  solution 
of  a  system  of  partial  differential  equations  (PDE) .  This  system  of  PDEs  is  in 
general  difficult  to  evaluate.  To  overcome  this  difficulty,  we  recognize  that  the 
PDE  system  solution  contributing  to  the  Hamiltonian  is  related  to  the  stochastic 
process  that  can  be  computationally  generated.  Consequently,  we  propose  to  use 
samples  of  the  stochastic  processes  to  evaluate  the  Hamiltonian  components.  In 
this  way,  our  approach  is  considerably  different  from  stochastic  optimal  control 
work  presented  in  [8] .  There  the  stochastic  processes  have  been  used  only  as  an 
analytical  tool  to  map  stochastic  process  to  be  controlled  into  the  finite  state 
space,  in  which  the  optimization  is  performed. 

The  benefit  of  using  solution  based  on  sampling,  i.e,  computational  statisti¬ 
cal  methods,  is  that  the  control  problems  in  robotics  could  be  solved  efficiently 
(in  real-time)  in  an  optimal  control  framework.  This  possibility  also  depends 
on  the  ability  to  implement  sampling  and  computations  with  samples  into  the 
processor  computing  the  control.  While  the  design  of  such  processors  is  a  sep¬ 
arate  issue,  ideally  it  should  be  driven  by  real  stochastic  processes  such  as  a 
noise  in  electronic  components.  This  concept  can  ultimately  provide  miniature 
computational  hardware  for  solving  complex  multi-robot  control  problems.  In 
this  respect,  the  ultimate  processor  will  be  a  system  in  which  computations  are 
based  on  laws  of  statistical  physics. 


In  Section  2  of  this  report  we  discuss  hybrid  systems  modeling  and  stochastic 
optimal  control  framework  applied  to  robotic  systems.  Section  3  introduces  a 
computational  method  for  computing  evolution  of  the  state  probability  density 
function  (PDF)  of  our  model.  In  this  section  we  compare  the  obtained  results 
with  the  solution  previously  obtained  from  the  PDF  system  describing  the  same 
evolution.  Section  4  describes  a  computational  method  for  evaluation  of  adjoint 
state  distributions.  Section  5  considers  a  special  case  of  the  process  self-adjoint 
systems.  Conclusion  and  suggestions  for  future  work  are  provided  in  the  last 
section. 
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2  Modeling  and  Control  Framework 


In  the  modeling  framework  we  consider,  the  state  of  an  individual  robot  at  time 
t  is  uniquely  defined  by  the  couple  {x{t),q{t)),  x  G  X,  X  G  M",  q  G  Q,Q  = 
{1,2,  While  in  the  discrete  state  (mode)  k  €  Q  the  continuous  state  of 

robot  obeys  differential  equation  x  =  fk{x,t).  We  also  assume  that  switching 
among  discrete  states,  say  from  the  state  k  G  Q  to  the  state  j  G  Q,  {k  ^  j) 
is  described  by  stochastic  transition  rates  \kj,  and  that  x{t)  is  a  continuous 
function  of  time.  In  other  words,  the  continuous  state  just  before  the  discrete 
state  transition  x(t~)  is  equal  to  the  state  x{t'^)  after  the  state  transition.  This 
very  general  model  of  individual  robot  is  illustrated  in  Fig.  1  and  the  modeling 
framework  we  are  applying  here  is  detailed  in  [10]. 

Recognizing  that  the  state  of  the  individual  robot  is  composed  of  discrete 
and  continuous  components,  the  state  probability  density  function  (PDF)  is 
a  vector  of  functions  p{x,t)  =  [pi{x,t)  p2{x,t) . . .  pN{x,t)]' .  Each  component 
Pi{x,t)  corresponds  to  the  discrete  state  i,  and  the  symbol  (’)  denotes  vector 


Figure  1:  Stochastic  hybrid  automaton  model  of  a  robot  in  the  probabilistic 
framework:  discrete  state  q;  continuous  vector  x  vector  field  fk,  k  G  Q  describes 
change  of  continuous  state  ;  stochastic  transition  rates  Xkj,  k,j  G  Q  describe 
mode  switching;  y  is  the  measurable  output,  if  the  full  continuous  state  of  the 
robot  is  measurable,  C  is  unity  matrix. 
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transpose.  The  state  PDF  satisfies 

/  Pi{x,t)dx  =  '^Pi{t)  =  1,  where  Pi{t)  =  /  pi{x,t)dx  (1) 
ieQ 

where  Pi{t)  is  the  probability  of  the  discrete  state  i  at  the  time  point  t.  Let  us 
define  the  vector  of  discrete  probabilities  P{t)  =  [Pi{t),  P2{t), . . .  PN{t)]' ,  then 
evolution  of  the  probability  vector  is  given  by: 

P{t)  =  Ft{t)P{t),  where  [Ft]ij  =  (2) 

with  matrix  Ft  defining  the  transition  rates  among  the  discrete  states.  In  gen¬ 
eral,  correspondence  between  the  matrix  Ft  members  [Ft]ij  and  the  transition 
rates  Xij  is  not  one-to-one.  Assuming  that  the  transition  rates  depend  on  a  vec¬ 
tor  u{t)  =  [ui{t)  U2(t)  . . .  UM{t)]'  of  variables  Ui,i  =  1,2...  M,  we  can  define 
the  transition  rate  matrix  as  a  function  of  the  vector  u{t),  i.e.,  Ft{t)  =  Fu{u{t)). 
Consequently,  the  vector  of  the  discrete  state  probabilities  obeys  [1]: 

P{t)  =  Fu{u)P{t)  (3) 


Moreover,  it  can  be  proven  [10]  that  the  state  PDF  obeys  the  following  system 
of  partial  differential  equations  (PDF) : 

=  F{u)p{x,  t)  =  {Fu{u{t))  +  Fo)p{x,  t)  (4) 


where  Fg  is  the  diagonal  linear  differential  operator.  When  the  operator  Fg  is 
applied  to  p{x,t),  it  results  in: 


\Fgp{x,  t)]ij 


V  •  {f^p^{x,t)),i  =  j 


(5) 
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Taking  into  account  that  the  state  PDF,  p{x,t),  evolution  depends  on  the 
vector  u{t)  we  can  formulate  the  optimal  control  problem  in  the  probability 
space  using  the  cost  function: 

J=  /  w' {x)p{x,T)dx  (6) 

J  X 

In  this  respect,  the  optimal  control  problem  is  the  optimization  problem: 

u*{t)=  max  J=  max  /  w{x)p{x,T)dx  (7) 

u(t)GUad  u{t)GUadJx 

Alternatively,  to  avoid  the  singular  control  problems  [10]  we  can  also  consider 
the  optimal  control  that  includes  the  term  penalizing  the  control: 

u*{t)=  max  J=  max  [  w{x)p{x,T)dx  +  e  [  u'{t)u'{t)  (8) 

u(t)^Ua.d  u{t)GUadJx  Jo 

Anyway,  the  solution  of  this  problem  is  a  sequence  of  the  optimal  control  u*(t), 
from  the  set  of  admissible  control  Uad,  such  that  the  cost  function  is  maximized. 
By  suitable  choice  of  the  weighting  function  w(x),  the  cost  function  can  be  used 
to  find  the  optimal  control  maximizing  probability  of  robot  presence  in  the 
desired  region  of  the  robots’  operating  space. 

The  optimal  control  maximizing  the  criterion  (6)  is  a  special  case  of  a  more 
general  optimal  control  problem  of  the  evolution  equation  [4].  Under  the  con¬ 
dition  that  the  operator  F(u)  is  bounded,  i.e.,  ||F'(M(t))||  <  oo  the  minimum 
principle  for  PDFs  can  be  applied  [4] .  According  to  the  minimum  principle,  the 
optimal  control  u*(t)  satisfies: 

u*(t)  =  arg  min  H {p{x ,  t) ,  u{t) ,  t)  (9) 

U&Uad 

In  other  words,  for  the  optimal  state  PDF  trajectory  p*{x,  t),  the  optimal  control 
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minimizes  the  Hamiltonian  at  each  time  point.  The  Hamiltonian  is: 


H{p{x,t),u,t)  =  {-K{x,t),F{u)p{x,t)) 


(10) 


where  brackets  (•,  •)  denote  the  scalar  product  of  function  vectors  defined  as: 


{p{x),q{x))  =  /  p'{x)q{x)dx 
Jx 


(11) 


The  function  vector  7r(x,  t)  is  the  so-called  adjoint  state  distribution  and  satisfies: 


dTr(x,  t) 
dt 

7r(x,r)  = 


—F' (u)tt{x,  t) 
—w(x) 


(12) 

(13) 


The  major  limitation  of  the  modeling  approach  presented  here  is  that  the 
system  of  PDEs,  derived  directly  from  the  model  shown  in  Fig.  1,  is  difficult 
to  solve  in  general  for  any  type  of  function  fi{x)  and  of  more  than  a  single 
dimension.  Nevertheless,  the  probabilistic  approach  we  use  is  of  fundamental 
importance  for  establishing  mathematically  tractable  relation  for  the  control  of 
spatio-temporal  distribution  of  large-size  robotic  populations.  Once  we  establish 
this  relation,  we  can  introduce  a  different  type  of  the  propagator  which  is  not 
necessarily  in  the  form  of  differential  equations.  Ultimately,  we  can  design  an 
optimal  control  taking  into  account  local  interactions  among  the  robots,  or 
between  the  robots  and  the  environment. 
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3  Stochastic  Sampling  Propagator 

Evolution  of  the  large-size  population  probability  density  function  p{x,  t)  is 
described  by  the  PDE  system  (see  Section  2,  Eq.4).  One  way  to  obtain  the 
evolution  p{x,  t)  is  to  solve  the  PDE  system  forward  in  time  starting  from  an 
initial  condition  p(x,  0)  =  p^{x).  This  numerical  solution  is  sensitive  to  function 
types  used  in  the  model,  initial  and  boundary  conditions  as  well  as,  in  finite 
element  approximation,  to  the  mesh  selection.  Since  numerical  errors  of  solu¬ 
tion  accumulate,  we  also  need  to  limit  the  terminal  time  T  of  the  solution.  If 
the  numerical  method  for  the  solution  is  not  carefully  designed  the  influence 
of  approximations  and  accumulated  errors  have  almost  random  effect  on  the 
solution. 

The  Hamiltonian  evaluation  is  instrumental  for  computing  the  optimal  con¬ 
trol  based  on  minimum  principle.  In  [9,  10]  we  formulated  the  algorithm  for 
the  optimal  control  in  ID  case.  However,  it  appears  that  sensitivity  of  numeri¬ 
cal  solutions  contributing  to  the  Hamiltonian  evaluation  is  a  limiting  factor  for 
solving  the  optimal  control  problem  in  more  than  one  dimension  and  with  larger 
number  of  discrete  states. 

Here  we  propose  an  approach  to  compute  the  evolution  p{x,  t)  based  on 
stochastic  trajectories  of  the  hybrid  state  {x,q)  evolution  resulting  from  the 
model  presented  in  Fig.  1,  Section  2.  To  account  for  the  fact  that  the  transition 
rates  can  change  in  time,  we  assume  that  the  control  is  a  piecewise  constant 
function  of  time  discretized  with  the  sample  time  AT.  The  basis  for  the  pro¬ 
posed  algorithm  is  the  Gillespie’s  stochastic  simulation  algorithm  [5]. 

To  generate  the  trajectory  of  (x,  q)  we  need  to  generate  the  initial  state 
(a;(0),  (7(0))  from  the  state  PDF  p{x,0)  =  p^{x).  Probability  Pi{t)  of  q{t)  =  i  is: 

Pi{t)  =  [  p,{x,t)dx  (14) 

Jx 
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Therefore,  the  random  variable  q{0)  =  i  should  be  generated  from  the  the 
discrete  state  probability  distribution  represented  by  the  vector  of  discrete  state 
probabilities  P(0)  =  [A(0)  -P2(0)  . . .  PAf(O)]'.  Symbolically,  we  will  represent  it 
as: 

9(0)=i-P(0)  (15) 

Once  the  initial  discrete  q{0)  state  is  defined,  the  continuous  variable  a;(0)  can 
be  generated  from  the  corresponding  pi{x,  0)  component  of  the  state  PDF  ,i.e., 
from  the  probability  V  of  x{t)  given  that  q(t)  =  i  and  t  =  0 

a;(0)  ~  Vixlqit)  =  i,t  =  0}  =  p^{x,  0)/Pi(0)  (16) 

Whenever  the  discrete  state  is  q{t)  =  i,  the  evolution  of  the  continuous  state 
X  obeys  x  =  fi{x).  Therefore,  generating  trajectory  {x{t),q{t))  reduces  to  the 
problem  of  generating  the  state  transitions  of  the  discrete  state  q{t).  Let  us 
assume  that  at  time  t  =  0,  0  G  [{k—l)AT,  kAT)  the  hybrid  state  is  {x{t),q{t)) 
the  time  instant  at  which  the  state  changes  tc  can  be  generated  based  on  the 
following  two  rules: 

(a)  tc  =  ts  +  tt,  tt  ~  3  ^  under  condition  that  p  <  kAT.  If 

the  condition  is  not  satisfied,  apply  rule  (b). 

(b)  tc  =  kAT  +  tt,  tt  ~  under  condition  that  tc  <  {k+  1)AT. 

If  the  condition  is  not  satisfied  increase  k  by  1.  Apply  rule  (b)  until  the 
condition  is  satisfied. 

These  two  rules  define  the  time  point  p  at  which  the  jump  from  discrete  state 
i  to  discrete  state  j  happens,  but  do  not  specify  the  variable  j.  The  state  j 
needs  to  be  sampled  from  the  discrete  state  probability  density  function,  i.e., 
from  the  probability  V  of  q{t^)  =  j  given  that  q{t)  =  i  provided  in  the  vector 
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of  the  discrete  probability  distribution  with  TV  —  1  elements: 


(17) 


N-l 


The  algorithm  just  described  can  be  used  to  generate  a  single  trajectory  for 
the  stochastic  model  shown  in  Fig.  1.  In  the  limit  of  large  number  of  samples, 
the  normalized  density  of  trajectory  points  will  correspond  to  the  solution  the 
PDF  system  given  by  Eq.4  in  Section  2.  In  this  respect,  the  stochastic  simulation 
is  a  computational  propagator  of  the  evolution  p{x,t)  and  we  can  denote  it  as: 


(18) 


To  illustrate  and  verify  the  algorithm  for  generating  stochastic  trajectories 
{x{t),q{t))  we  use  the  following  ID  and  2D  examples. 

3.1  ID  example 

The  stochastic  model  presented  in  Fig.  2b  illustrates  the  state  PDF  evolution 
of  the  large-size  robot  population  along  one  dimension,  Fig.  2a  in  which  ui, 
U2  and  U3  correspond  to  stochastic  rates  of  commands:  move-left,  move-right 


q=l 


q=2 


a) 


Figure  2:  ID  example 
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and  stop.  In  this  example  ki  =  —0.5  and  ^2  =  0.25.  The  control  u{t)  = 
U2{t)  usit)]  is  computed  as  the  optimal  control  based  on  the  Minimum- 
principle  and  Hamiltonian  presented  in  the  previous  section. 

The  cost  function  is: 


J  = 


lx 


w'{x)p{x,t)dx  +  e  /  Ui{t)  +  ul{t)  +  u^{t)dt 


(19) 


where  £  =  10  the  weighting  w{x)  =  [0  0  W3(a;)]'  and  initial  condition  p{x,  0)  = 
[0  0  Psix,  0)]'  are  defined  by: 


wsix)  = 


1.25  <x<  2.25 
0,  otherwise 


(20) 


P3(a:,0)  = 


Vo. 0277 


rp(- 


(a:-2.5)" 

0.02 


),  2  <  a;  <  3 


0,  otherwise 


(21) 


The  optimal  control  sequence  u*{t)  =  [^^(t)  U2it)  u’^it)]  in  the  time  interval 
0  <  t  <  3  is  defined  by: 


ulit)  = 


2,  ti<t<t2 
0,  elsewhere 


ulit)  =  0,  ulit)  = 


2,  ti  <  t  <  t2 
0,  elsewhere 


(22) 


The  evolution  of  the  state  PDF  for  this  system  under  the  control  u*(t)  is  pre¬ 
sented  in  Fig.  3.  We  present  only  pi(x,  t)  and  psix,  t)  because  under  this  control 
P2ix,t)  =  0,  Vt. 

For  the  illustration  we  generated  10  stochastic  trajectories  of  continuous 
variable  x  (see  Fig.  4)  under  the  control  u*it).  Evolution  of  the  discrete  state 
q  can  be  observed  from  the  trend  in  x.  When  x  decreases  the  discrete  state  is 
1,  and  when  it  stays  constant  the  state  is  g  =  3.  It  is  worth  mentioning  that 
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0  0.5  1  1.5  2  2.5  3  3.5  4  0  0.5  1  1.5  2  2.5  3  3.5  4 

Position  X  [km]  Position  x  [km] 

Figure  3:  The  finite  element  solution  of  the  state  PDF  evolution  for  ID  example 
under  the  optimal  control  u*{t) 


Figure  4:  Random  set  of  10  trajectories  resulting  from  the  stochastic  simulation 
under  the  optimal  control  u*{t) 


among  these  10  trajectories  there  is  one  for  which  x{t)  is  constant.  The  small 
pick  at  around  the  point  2.5  in  the  right  panel  of  the  Fig. 3  at  t  =  3,  confirms 
that  the  probability  of  such  trajectories  is  non-zero,  but  it  is  small. 

To  obtain  the  state  PDF  p{x,t),  i.e.,  its  components  pi{x,t)  at  specific  time 
point  t  we  need  to  collect  points  x{t)  and  estimate  components  pi{x,t).  It 
is  obvious  that  10  trajectories  cannot  provide  a  good  estimate  of  p{x,t)  For 
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this  reason,  we  generated  10®  trajectories  and  computed  histogram  probability 
density  function  estimate.  That  means  that  we  discretized  the  x  axis  into  inter¬ 
vals  of  the  length  Acc  and  counted  how  many  points  fell  into  a  specific  region. 
Finally,  we  normalized  the  histogram  so  that  the  estimated  p{x,  t)  satisfied  con¬ 
dition  1,  given  in  Section  2.  Results  are  presented  in  Fig.  5.  As  expected  the 
match  between  the  numerical  PDF  system  solution  and  the  result  obtained  from 
stochastic  trajectories  is  exact.  There  are  only  negligible  discrepancies  due  to 
data  sampling  from  finite  number  of  trajectories. 


0  0.5  1  1.5  2  2.5  3  3.5  4  0  0.5  1  1.5  2  2.5  3  3.5  4 


Position  X  [km]  Position  x  [km] 

Figure  5:  The  stochastic  simulation  solution  of  the  state  PDF  evolution  for  ID 
example  under  the  optimal  control  u*{t) 

3.2  2D  example 

The  benefit  of  applying  stochastic  sampling  algorithm  and  the  reason  for  intro¬ 
ducing  the  propagator  is  its  ability  to  propagate  multi-dimensional  probability 
density  function.  To  illustrate  that,  we  apply  stochastic  sampling  for  computing 
distribution  of  the  model  presented  in  Fig.  6.  The  model  is  based  on  a  scenario 
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in  which  a  population  of  small  robots  is  initially  concentrated  on  a  location  in 
the  operating  region. 

The  robots  are  controlled  by  the  stochastic  signals  produced  by  aerial  robots 
(Fig.  6a).  Each  robot  in  the  population  moves  in  the  direction  of  the  active 
signal  source.  Under  the  assumptions  that  aerial  robots  are  far  away  from  the 
population  and  that  the  robot  velocity  is  constant  (v  =  1),  the  robotic  motion 
model  is  given  by  the  stochastic  model  shown  in  Fig.  6c.  In  this  example  no 
direct  transitions  between  states  1  and  3  exist.  The  model  and  scenario  are 
detailed  in  [10,  11] 


2D  example 

We  take  the  case  in  which  A12  =  0.5,  A21  =  0.1,  A23  =  0.9  and  A32  =  0.1 
are  constant.  The  trajectories  are  generated  based  on  the  proposed  algorithm 
and  sampled  at  time  points  t=0,  0.39,  0.79,  1.18,  1.57,  1.96.  The  samples  are 
presented  in  Fig.  7  and  the  color  of  dots  denotes  the  time  point  of  the  sample. 
To  produce  this  result,  we  used  only  100  samples  so  that  difference  in  the  density 
of  points  could  be  observed. 

The  transition  rates  of  this  example  are  set  without  consideration  of  optimal 
control.  Ultimately  we  should  be  able  to  control  these  2D  distributions  similar 
to  the  ones  in  the  ID  example  presented  in  the  previous  section.  For  comparison, 
we  also  present  in  Fig.  8  the  state  PDF  evolution  resulting  from  the  solution  of 
the  corresponding  PDE  system.  This  evolution  can  help  us  recognize  possible 
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Figure  7:  The  stochastic  simulation  solution  of  the  state  PDF  evolution  for  2D 
example  (100  points).  The  transition  rates  are  A12  =  0.5,  A21  =  0.1,  A23  =  0.9 
and  A32  =  0.1.  The  time  points  are  t=0,  0.39,  0.79,  1.18,  1.57,  1.96  and  are 
represented  by  different  colors  from  pink  to  black;  dimensions  along  x  and  y 
axis  are  in  km. 


but  improbable  stochastic  trajectory  realizations  presented  in  Fig.  7.  We  can 
also  note  the  polygon  like  contour  plots  resulting  from  the  imprecision  of  the 
numerical  solution  of  the  PDF  system. 


Figure  8:  Numerical  PDF  system  solution 
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4  Stochastic  Sampling  for  Adjoint  State  Distri¬ 


bution 

The  Hamiltonian,  corresponding  to  the  optimal  control  problem  Eq.7,  Section 
2,  is  defined  by: 


H{t)  =  /  TT{x,t)F{u)p{x,t)dx  =  Ep{T:{x,t)F{u)}  (23) 

Jx 


and  can  be  seen  as  an  expected  value  of  the  adjoint  state  (or  co-state)  with 
respect  to  the  the  state  PDF.  Assuming  that  the  adjoint  state  evolution  TT{x,t) 
is  known,  evaluation  of  the  expected  value  is  easy  having  in  mind  that  we  know 
how  to  generate  realizations  of  the  hybrid  state  stochastic  trajectory  {x{t),q{t)). 
One  way  to  obtain  the  adjoint  state  distribution  evolution  tt{x,  t)  is  to  solve  the 
PDF  system  Eq.  12,  Section  2,  backward  in  time  starting  from  the  initial 
condition  tt(x,T)  =  —w{x).  However,  we  should  notice  that  the  complexity  of 
solving  the  adjont  PDF  system  is  as  complex  as  solving  the  state  PDF  evolution 
p{x,  t).  Therefore,  it  is  natural  to  consider  an  opportunity  to  solve  adjoint  state 
evolution  7r(x,  t)  using  stochastic  sampling  similar  to  the  one  presented  in  the 
previous  section. 

The  adjoint  set  of  equation  is  defined  by: 


Tr(x,T)  =  —w(x) 


(24) 

(25) 


Let  us  assume  that  the  7r(x,t)  as  well  as  w{x)  are  from  the  state  PDF  class  of 
functions,  i.e.,  they  satisfy  the  following  property  {Property  1): 


ieQ 


(26) 
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(27) 


where  P^{t)  is  the  probability  of  discrete  state  i  at  time  t  assuming  that  the 
Tr{x,t)  is  the  state  PDF  of  appropriately  defined  stochastic  process.  If  the 
assumption  is  valid  than  the  discrete  state  probability  evolution  backward  in 
time  is  defined  as: 

p-(r-r)=F'(u)P"(T-T)  (28) 

but  here  we  find  the  contradiction.  The  transpose,  FJ),  of  the  transition  rate 
matrix  F„  does  not  warranty  in  general  Property  1.  Therefore,  the  adjoint  PDF 
system  does  not  define  a  conservation  law  in  the  probability  space.  This  is  a 
major  obstacle  for  the  straight-forward  design  of  an  algorithm  that  will  evaluate 
adjoint  state  distribution  evolution  based  on  the  stochastic  process  similar  to 
the  one  for  evaluating  the  state  PDF.  Consequently,  an  alternative  approach 
needs  to  be  considered. 

The  algorithm  for  computing  solution  of  adjoint  equation  has  been  recently 
published  in  [14].  The  algorithm  is  based  on  an  assumption  that  the  PDE 
system  describing  the  adjoint  state  distribution  evolution 

^  = -F'{u)TT{x,t)  (29) 

can  be  discretized  in  space  and  time,  so  that  we  obtain  a  system  of  linear 
equations  defining  the  solution  of  Eq.  29: 

Alt  =  h  (30) 

where  tt  is  a  vector  of  values  irAmAX,  kAT),  i  G  Q,  AT  is  the  discretization 
sample  time  and  m  corresponds  to  the  mth  element  of  the  space  discretization. 
Obviously  the  matrix  A  is  a  square  matrix  of  dimension  (|M|-|(5|-Fr)x(|M|- 
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IQI  •  K),  where  \M\,  \K\  and  \Q\  are  the  number  of  discrete  element  of  space  X, 
the  number  of  the  time  steps  and  the  number  of  the  discrete  states,  respectively. 
Vector  b  is  equal  to  0  except  for  the  elements  corresponding  to  the  time  T  where 
b  element  values  correspond  to  —w(x). 

Equation  29  is  an  algebraic  equivalent  of  Eq.  30.  Based  on  this  algebraic 
interpretation,  tt  can  be  obtained  as: 

TT  =  A~^b  (31) 

Instead  of  computing  the  inversion  algebraically,  the  authors  of  [14]  proposed 
an  algorithm  in  which  Eq.  30  is  rewritten  in  the  following  form: 

TT  =  Att  +  b,  A  =  I  —  PA,  and  b  =  Pb  (32) 

with  P  a  being  block-diagonal  preconditioner  matrix.  Based  on  this  form  the 
solution  for  tt  can  be  expanded  in  a  Neumann  series: 

TT  =  b  + Ab+ A^bP  A^b+ . . .  (33) 

This  series  converges  only  if  the  spectral  radius  of  A  is  less  than  1.  In  this  case 
we  can  use  the  Monte  Carlo  method  to  sample  this  infinite  series  by  a  Markovian 
random  walk. 

The  algorithm  we  present  here  is  adopted  from  [14]  and  the  specific  values 
are  adjusted  to  the  special  case  of  dynamical  system  we  are  dealing  with.  The 
single  realization  on  the  random  walk  we  denote  by  [p\.  The  value  lT[p],  or  D[p] 
denotes  the  variable  W,  or  D  related  to  the  single  realization  of  the  random 
walk.  Variable  ak\p]  denotes  the  state  of  the  random  walk  at  the  step  k  from 
the  realization  of  the  random  walk  [p] .  Introducing  this  notation,  the  algorithm 
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to  compute  tt  is  as  follows: 


1.  Initialize  W[p]  =  1,  D[p]  =  0  and  set  t  =  1 

2.  Compute  the  matrix  A,  and  based  on  it,  define  probabilities  of  state  tran¬ 
sitions: 

P{hj)  =  ^"^'4  (34) 

3.  For  each  random  walk  \p]  that  is  at  time  t  choose  the  next  state  based  on 
random  probabilities  p{i,j).  If  the  state  ak+i[p]  is  not  the  final  exit  state, 
update  W[p]  and  D[p]  as  follows: 

W[p]  =  IF[p]wa,b]afc+ib];  D[p]  =  D[p]  +  (35) 

If  ak+i[p]  is  the  final  exit  state,  the  random  walk  is  absorbed  and  we  freeze 
D[p\. 

4.  Repeat  state  3  until  all  random  walks  at  time  step  t  are  either  absorbed 
or  left  the  time  step.  If  t  <  m,  then  let  t  =  t  -I-  I  and  go  to  step  2. 

5.  After  being  done  with  the  last  time  step  m,  all  random  walks  are  ab¬ 
sorbed.  Compute  the  sample  mean  of  the  estimator  which 

is  approximation  of  tt. 

Evaluation  of  the  state  PDF  as  well  as  adjoint  state  distribution  based  on 
stochastic  simulation  provides  an  opportunity  for  the  time  efficient  computing  of 
the  optimal  control.  However,  computations  need  to  be  carefully  analyzed  since 
the  adjoint  state  evolution  is  not  a  conservation  law  and  it  can  result  in  com¬ 
putational  instability.  Hopefully,  a  sacrifice  of  precision  of  the  adjoint  state  dis¬ 
tribution  computations  would  not  influence  much  the  Hamiltonian  evaluation. 
Based  on  our  experience,  the  state  PDF  evaluation  presented  in  the  previous 
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section  is  much  faster  than  the  evaluation  based  on  solving  partial  differential 
equations.  Similar  tests  should  be  done  with  the  computation  of  the  adjoint 
state  distributions  presented  in  this  section.  With  the  design  of  dedicated  hard¬ 
ware,  we  believe  that  all  this  difficulty  can  be  avoided  and  ultimately  complete 
stochastic  sampling  decision  making  can  be  done  on-board  and  in  real-time. 

5  Process  Self-adjoint  Systems 

In  the  previous  section  we  explained  that  in  general  the  adjoint  state  equation 
is  not  a  conservation  law.  Consequently  the  stochastic  sampling  computation 
of  adjoint  state  distributions  is  a  more  difficult  job  than  estimation  of  the  state 
PDFs.  In  this  section  we  consider  a  special  type  of  dynamical  systems  in  which 
the  adjoint  state  evolution  is  a  conservation  law,  and  we  call  this  type  of  systems 
process  self-adjoint. 

Let  us  assume  that  that  the  state  pdf  evolution  is  described  by: 

=  F'{u)p{x,t)  (36) 

p{x,0)  =  p°{x)  (37) 

and  the  adjoint  state  distribution  evolution  is  given  by: 

=  -F'{u)Tr{x,t)  (38) 

tt{x,T)  =  —w{x)  (39) 

where  the  operator  F{u)  =  {Fu{u{t))  +  Fg)  and  the  operator  F'{u)  is: 

F'{u)  =  {F^{u{t))  +  Fg)'  =  {F'^{u{t))  -  Fg)  (40) 
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In  other  words,  the  adjoint  state  distribution  evolution  is  rewritten  as: 


-(Kiuit))  -  Fd)TT{x,  t)  (41) 

—w{x)  (42) 

If  we  rewrite  this  equation  backward  in  time  we  obtain: 

{Kiuit))  -  Fd)TT{x,  T-t)  (43) 

—w{x)  (44) 

Finally,  introducing  the  substitution  ^(t)  =  —Tr(x,  T  —  t)  we  can  obtain  the 
backward  solution  of  the  adjoint  state  evolution  in  the  form  of  a  forward  time 
PDF  system: 


d'K{x,  T  —  t) 
dr 

7r(a:,  r  =  0)  = 


dTT{x,  t) 
dt 

tt{x,T)  = 


d(l){x,T) 

dr 

(l}{x,0) 


(Kiuit))  -  Fd)4>{x,T) 
w{x) 


(45) 

(46) 


Let  us  consider  the  special  case  in  which  the  discrete  state  transition  matrix 
satisfies  Fl^{u{t))  =  F„(u(t)).  Then  we  can  conclude  that  the  adjoint  state  evo¬ 
lution  Eq.  5  (i.e.,  its  equivalent  backward  time  evolution  Eq.  5)  is  conservation 
law.  Consequently,  if  w{x)  is  a  function  of  the  same  type  as  the  state  PDF, 
then  the  adjoint  state  distribution  tt,  i.e.,  its  equivalent  0  can  be  considered  the 
co-state  probability  density  function  of  some  stochastic  process  which  we  will 
explain  briefly. 

The  condition  F!^{u{t))  =  Fu{u{t))  implies  that  the  transition  rate  from  the 
discrete  state  i  to  the  discrete  state  j  is  the  same  in  the  opposite  direction,  i.e., 
from  the  state  j  to  the  state  i.  Interestingly  enough,  under  assumption  that 
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=  Fu{u{t))  the  robotic  system  can  be  considered  process  self-adjoint. 
In  other  words,  the  adjoint  state  evolution  corresponds  to  a  stochastic  process. 
For  example,  let  us  consider  the  vehicle  of  which  dynamics  can  be  represented  by 
three  discrete  and  three  continuous  states  (see  Fig.  9a  and  vehicle  illustration 
in  Fig.  10).  The  first  discrete  state  corresponds  to  the  vehicle  moving  without 


Figure  9:  a)  The  stochastic  process  describing  the  evolution  of  the  state  PDF, 
b)  The  stochastic  process  describing  the  evolution  of  the  adjoint  state  PDF 
backward  in  time  (j) 

change  of  direction.  The  second  and  the  third  state  correspond  to  the  increase 
and  decrease  of  the  angle  of  the  moving  direction,  respectively.  To  obtain  a 
stochastic  process  describing  evolution  of  the  adjoint  state  distributions  we  need 
each  dynamics  of  the  discrete  state  i,  fi  substitute  with  —fi  and  keep  in  mind 
that  this  stochastic  process  describes  the  adjoint  state  backward  in  time.  The 
resulting  adjoint  process  is  depicted  in  Fig.  9b.  The  generation  of  this  stochastic 
process  can  be  performed  using  the  algorithm  presented  in  Section  3. 

To  summarize,  evaluation  of  the  state  p{x,t)  and  co-state  Tr{x,t)  PDFs  for 
this  case  will  be  based  on  stochastic  samples  from  trajectories  defined  by  models 
presented  in  Fig.  9.  The  trajectories  need  to  be  generated  and  collected  in  3D 
since  we  deal  with  the  three  continuous  state  variables.  The  2D  projection  of 
these  3D  trajectories  is  illustrated  in  Fig.  10.  The  distribution  of  initial  points 
of  trajectories  generated  forward  in  time  relates  to  p^{x)  in  the  region  A  (see 
Fig.  10)  while  the  initial  distribution  of  trajectories  generated  backward  in 
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Figure  10:  a)  The  vehicle  and  continuous  state  variables  PDF,  b)  The  forward 
and  backward  paths  starting  from  the  region  A  and  B,  respectively 

time  relates  to  w{x)  in  the  region  B.  In  this  respect  the  Hamiltonian  Eq.  23 
measures  the  match  between  forward  and  backward  trajectories  and  its  value 
can  be  used  for  control  improvement,  i.e.,  to  provide  a  better  match.  Although 
we  discuss  the  special  case,  this  is  an  elegant  example  providing  rationale  behind 
the  ant-optimization  [3]  algorithms. 

6  Conclusions  and  Future  work 

In  this  research  project  we  study  a  theoretical  possibility  of  solving  the  opti¬ 
mal  control  problems  based  on  stochastic  sampling.  We  came  to  the  conclusion 
that  the  main  difficulty  in  this  approach  is  evaluation  of  the  adjoint  state  dis¬ 
tributions.  While  proposed  Monte  Carlo  algorithm  is  numerically  efficient,  its 
convergence  is  not  guaranteed.  An  alternative  approach  to  avoid  the  difficulty 
is  to  consider  a  class  of  robotic  systems  for  which  the  adjoint  state  distribution 
evolution  corresponds  to  the  stochastic  process.  We  call  such  systems  process 
self-adjoint.  In  this  case  both  state  PDF  and  co-state  PDF  contributing  to 
the  Hamiltonian  can  be  evaluated  based  on  stochastic  samples.  Based  on  the 
process  self-adjoint  systems  we  can  establish  a  relation  between  the  stochastic 
sampling  stochastic  optimal  control  and  ant-optimization  algorithms  [3]. 

Future  work  along  this  approach  will  be  focused  on  examples  based  on  algo- 
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rithms  presented  here.  We  are  also  considering  establishing  relation  between  the 
stochastic  sampling  methods  and  Hamilton- Jacobi-Bellman  equation  in  order  to 
solve  the  stochastic  optimal  control  problems  in  a  closed  loop. 

The  main  idea  behind  our  work  is  to  solve  stochastic  optimal  control  prob¬ 
lems  based  on  stochastic  sampling.  We  believe  that  in  this  way  we  can  employ 
modern  statistical  methods  and  control  theory  to  solve  complex  stochastic  op¬ 
timal  control  problems  for  multi-robot  systems  in  real-time. 
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Abstract 

Motility  is  an  important  property  of  immune  system  cells.  It  provides  cells 
with  the  ability  to  perform  their  function  not  only  at  the  right  time,  but  also 
at  the  right  place.  In  this  paper,  we  introduce  the  problem  of  modeling  and 
estimating  an  effective  force  field  directing  cell  movement  by  the  analysis  of 
intravital  video  microscopy.  A  computational  approach  is  proposed  for  solving 
this  problem  without  dealing  with  a  parameterized  spatial  model  of  the  field 
in  order  to  avoid  potential  errors  due  to  inaccurate  spatial  model  assumptions. 
We  consider  the  dynamics  of  cells  similar  to  the  dynamics  of  distributed  agents 
typically  used  in  the  field  of  swarm  robotics.  The  method  utilizes  a  fixed- 
interval  Kalman  filter  based  smoother.  Its  application  results  in  a  map  giving 
the  intensity  and  direction  of  the  effective  force  field.  The  results  show  that  real¬ 
time  video  images  are  a  source  of  data,  enabling  us  to  visualize  intriguing  spatio- 
temporal  phenomena  inside  immune  system  organs.  The  proposed  approach  can 
fill  the  existing  gap  between  contemporary  technology  and  quantitative  data 
analyses  present  in  the  field  of  biosystems. 
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FORCE  LOCALIZATION  AND  MAPPING  FROM  INTRAVITAL 

VIDEO  MICROSCOPY  4 


Original  microscopy  video 


Dejan  Milutinovic,  Devendra  P.  Garg 
Robotics  and  Manufacturing  Automation  Laboratory,  Duke  University 

•  Video  microscopy  provides  information  about  immune  system  cells  residing  in  tissues  of  living  organisms. 

•  Force  Localization  And  Mapping  (FLAM)  is  the  problem  of  estimating  an  effective  force  field 
influencing  cell  motility  from  intravital  video  microscopy. 

•  For  solving  the  FLAM,  we  propose  a  computational  approach  without  dealing  with  a  parameterized 
spatial  model  of  the  field.  In  this  way,  we  avoid  potential  errors  due  to  inaccurate  spatial  model  assumptions. 

•  The  estimation  is  based  on  the  individual  agent  motility  model. 

•  No  alternative  microscopy  method  exists  for  measuring  these  forces 

•  The  computational  method  does  not  interfere  with  the  physiological  condition  inside  organs. 
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Figure  from  D.  Milutinovic.  P.  Lima,  "Cells  and  Robots",  Springer,  2007 


Force  profile: 

horizontal  component  constant  Mn  this  example,  we  can  easily  identify 


vertical  component  sinusoidal 


the  horizontal  force  component,  but  the 
analysis  does  not  provide  an  insight  into 
V^the  vertical  force  component. _ 


Single  cell  track  estimation  performance 
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(s:  the  Laplace  domain  variable) 


Data  are  known;  this  is  a  smoothing, 
not  a  filtering  problem 


Algorithm 


1 .  Estimate  the  force  as  if  a  constant  for  each  trajectory  I 

(  Smoother =  0) 

2.  “Average"  all  estimates  to  obtain  the  constant  force  term 


3.  Estimate  the  force  along  each  i‘^  trajectory  initially  assuming 
that  it  has  value  of  the  constant  force  obtained  in  Step  2  above 

( Smoother *0) 

4.  “Average”  all  estimates  that  belong  to  the  same  grid  cell  ••  ► 


Trajectory  lengths  and  estimation  fusion 


/rhe  cell  tracks  are  short  in  comparison  to  the  number 
of  data  points  necessary  for  reliable  force  estimation; 
We  fuse  estimations  corresponding  to  the  positions  of 
the  same  grid  eiement. 


Conclusions 
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♦  An  optimal  smoother  applied  to  force  estimations  as  sensed  by  individual  cells. 

♦  Estimations  averaged  over  a  grid  of  rectangular  regions  taking  into  account  the 
estimation  covariances. 

♦  The  force  field  is  estimated  and  visualized  by  the  average  forces  assigned  to  each 
rectangular  region.  No  analytical  constraints  on  the  force  field. 

♦  To  apply  the  method  to  the  analysis  of  imaging  data  routinely,  a  rapid  adjustment  of 
the  estimator  parameters  is  important. 

♦  Our  results  show  that  real-time  video  images  are  a  source  of  data,  enabling  us  force 
visualization,  which  can  be  used  for  studying  intriguing  spatio-temporal  phenomena 
inside  immune  system  organs. 

♦  The  problem  formulation  and  the  solution  fill  the  existing  gap  between  contemporary 
technology  and  quantitative  data  analyses  present  in  the  field. 


